OCS in small para-hydrogen clusters: energetics and structure with A^=l-8 complexed 

hydrogen molecules 
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I. INTRODUCTION 



F. Paesani, R. E. Zillich and K. B. Whaley 
Department of Chemistry and Pitzer Center for Theoretical Chemistry, University of California, Berkeley, CA 94720 

We determine the structure and energetics of complexes of the hnear OCS molecule with small 
numbers of para-hydrogen molecules, A''=l-8, using zero temperature quantum Monte Carlo meth- 
^ ' ods. Ground state calculations are carried out with importance-sampled rigid body diffusion Monte 

' Carlo (IS-RBDMC) and excited state calculations with the projection operator imaginary time spec- 

, tral evolution (POITSE) methodology. The ground states are found to be highly structured, with a 

gradual build up of two axial rings as A'' increases to 8. Analysis of the azimuthal density correla- 
tions around the OCS molecule shows that these rings are quite delocalized for small A'' values, but 
become strongly localized for N > 5 . Excited state calculations are made for a range of total cluster 
angular momentum values and the rotational energy levels fitted to obtain effective rotational and 
distortion constants of the complexed OCS molecule as a function of cluster size A'^. Detailed anal- 
ysis of these spectroscopic constants indicates that the complexes of OCS with para-hydrogen have 
an unusually rich variation in dynamical behavior, with sizes A^=l-2 showing near rigid behavior, 
C/3 , sizes A'^=3-4 showing extremely floppy behavior, and the larger sizes A'^=5-8 showing more rigid 

^ ' behavior again. The large values of the distortion constant D obtained for A'^=3-4 are rationalized 

^ I in terms of the coupling between the OCS rotations and the "breathing" mode of the first, partially 

^ . filled ring of para-hydrogen molecules. 

> 

O 

Complexes of small molecules with variable number of para- hydrogen (P-H2) molecules pose a challenging and 
^►-j fascinating arena for the manifestation of nuclear quantum effects on structure and spectroscopy. Like their better 
known analogous helium clusters, para-hydrogen clusters are expected to be strongly influenced by the quantum 
nature of the P-H2 molecules. The mass of P-H2 being one half that of ''He actually suggests that larger zero point 
effects might be expected. However, the greater binding of P-H2 to itself and to other species means that the greater 
' delocalization tendency of P-H2 competes with greater localizing potential energy terms. The competition between 
^ these two strong effects is what causes both calculation and analysis of para-hydrogen clusters and complexes to 
' be considerably more difficult than that for analogous helium systems. Like helium, para-hydrogen is a boson, but 
C , because of the greater binding energy, p-B.2 solidifies in the bulk before cooling to a low enough temperature for 
macroscopic manifestation of the boson permutation exchange symmetry to occur with a superfluid phase. However, 
in finite clusters vestiges of superfluidity may occur, particularly when the para-hydrogen packing is expanded or 

■ otherwise constrained by the reduced dimensionality and/or the 15-H2 molecules are bound to foreign species. 

■ The linear OCS molecule has played a key role in developing our understanding of the fundamental forces controlling 
the properties of P-H2 complexes. Experimental studies with OCS complexed by variable numbers of para-hydrogen 

^ (N < 16) and embedded in helium droplets have shown a very rich spectroscopic behavior, with both vibrational shifts 
' ^ and rotational fine structure providing indirect information on the symmetry and distribution of the P-H2 around the 
OCS molecule [1-3]. In helium, the rotational spectra of the para- hydrogen complexes are well fit by asymmetric 
^ tops for A^=l-4, while larger clusters appear adequately described by a symmetric top Hamiltonian [2]. For certain 
Qh' sizes, iV=5, 6, 11, 14-16, no Q-branch is seen in the spectra. For the larger sizes, A^=14-16, this disappearance is 
' temperature dependent, which led to the suggestion that it might be a manifestation of superfluidity in the para- 
hydrogen component. Path integral Monte Carlo calculations for a complete solvation shell (A^=:17) have shown that 
a transition to an anisotropic superfluid state is indeed found at temperatures below T= 0.3 K, and that this can 
account for the disappearance of the Q-branch for > 11 [4]. However, for the smaller sizes, the spectral anomaly 
is independent of temperature. In the absence of detailed knowledge of the structure and rigidity of the complex, the 
microscopic origin of this anomaly is less clear, although the permutation symmetry can also be expected to play a 
role. 

Much less is known experimentally about the complexes of OCS with para-hydrogen in the absence of a solvating 
helium environment. The 0CS(p-H2) dimes has been characterized by high resolution spectroscopy [5], and its 
strcuture has been shown to be approximate with that of an asymmetric top rotor. The structure of this dimer is 
very similar in the gas phase and in a helium droplet [6,7]. For complexes with larger numbers of para- hydrogen 
molecules, no spectroscopic data have been published. We expect that OCS complexed by 1 < iV < 17 para-hydrogen 
molecules will show a very rich variation in structure, energetics, and spectroscopic properties, as a consequence of the 
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strongly modulated solvating hydrogen density along the axis of the OCS molcciile seen in path integral calculations 
for TV =17 [4]. These modulations are considerably stronger than those found in the analogous complexes of OCS 
with helium [8-10] because of the stronger (p-H2)-0CS and (p-H2)-(p-H2) interactions. Consequently we can expect 
a considerably more complex size dependence of the spectroscopic properties for the complexes with p-H2 as the first 
solvation shell is filled. The small complexes of OCS with helium have recently been shown to undergo a transition 
from near rigid molecular complex behavior to true quantum (or "superfluid") solvation characterized by permutation 
exchanges along the molecular axis, as the number of helium atoms increases from A''=l to A'^=20 [10]. Given the recent 
demonstration of anisotropic superfluidity in OCS(p-H2)i7 [4], we expect an analogous transition for complexes with 
molecular hydrogen. However, the considerably greater competition between quantum (kinetic) delocalization and 
(potential) localization for para-hydrogen renders even qualitative prediction of the spectroscopic behavior impossible 
without detailed microscopic calculations. Furthermore, for the intermediate sizes 1 < A'' < 17, the properties of gas 
phase OCS(p-H2)i7 complexes and their analogous embedded in helium droplets can be expected to show maximal 
differences, since this is the size regime in which subtle changes in para-hydrogen and helium density distributions 
can give rise to significant energy differences. 

In this paper we present results of ground and excited state calculations for OCS complexed with N=l-8 para- 
hydrogen molecules, as part of a larger systematic study of the complexes with up to a complete solvation shell [11]. 
We employ two kinds of zero temperature quantum Monte Carlo methods, the diffusion Monte Carlo (DMC) and the 
projection operator imaginary time spectral evolution (POITSE) approach. These microscopic calculations provide 
energy levels and rotational and distortion constants, in addition to energetics and ground state structures. We find 
a very interesting variation in structure for this series of small complexes, showing considerably more inhomogencous 
behavior than the analogous series of complexes of OCS with ^He. We determine the structure both along the 
molecular axis and around it, with analysis of para- hydrogen pair correlation functions and find that this allows us to 
provide a good description of the extent of localization in the ground state. The excited state energies arc extracted 
from a maximum entropy analysis of the imaginary time correlation functions derived from the POITSE approach. 
We find that a simple multi-exponential fit of these correlations functions provides similar results, but with larger 
uncertainties. Analysis of the fitted spectral constants resulting from the excited state energies reveals a gradual 
change of symmetry of the 0CS(p-H2)jv complexes as A'' increases, which is accompanied by a monotonic decrease in 
the effective rotational constant over this size range. While overall the analysis of spectral constants indicates more 
rigid structures for these OCS complexes with para-hydrogen than those with helium, the fitted spectral constants 
nevertheless show a markedly strong appearance of floppy, non-rigid structures for both the N=3 and N=4: complexes. 
We discuss these structures and predicted spectroscopic constants in relation to available experimental data and in 
the context of the analogous OCS('*He)jv complexes [10]. 

II. THEORETICAL APPROACH 
A. Interaction Potentials 

The total interaction potential for the 0CS{p-B.2)n clusters is obtained as a sum of all pairwise contributions, 

neglecting the effects of three-body forces among the para-hydrogen molecules, while correctly treating the many- 
body interaction between the linear OCS impurity molecule and each P-H2 molecule one at a time. Because of its 
nuclear spin symmetry (1=0), para- hydrogen can only access even values of its angular momentum, i.e., j/f=0, 2, 4, 
etc. Furthermore, given the large spacing between the rotational levels, jn can generally be assumed to be a good 
quantum number in the weakly bound van der Waals complexes [12,13,5]. These considerations, together with the 
fact that our calculations arc carried out at a temperature of K, allow us to consider each P-II2 molecule as being in 
its jH=0 rotational state. Therefore, in the present study, we treat the para- hydrogen molecules as spherical particles. 
The total interaction potential can then be expressed as 

N 

v{R) = ^y°^^-(f-^=)(i?i,i?i) + ^y(p-f^=)-(f-f^=)(iii,) (1) 

i=l i<j 

where R is a generic vector that defines the coordinates of all the particles. Therefore, Rij is the distance between 
P-H2 molecules i and 7, and "di are the Jacobi coordinates of the ith para- hydrogen molecule in the center of mass 
frame of the OCS (see Fig. 1). 

For the (p-H2)-(p-Il2) interaction we use the spherical part of the empirical potential proposed by Buck et al. [14]. 
The 0CS-(p-H2) interaction is obtained after integration over the P-H2 angular variables (i?',0') of a 4-dimensional 
ah initio potential energy surface recently calculated by Higgins et al. using fourth-order MoUer-Plesset perturbation 
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theory [15]. A contour plot of the averaged potential is shown in Fig. 2a. The global minimum of -144.5 cm~^ is 
located at R ~ 3.35 A and i9 = 105°. Two other local minima corresponding to the two coUinear geometries are also 
found. The first, at i9 = 0° and R = 4.52 A has a well depth of -91.4 cm~^ The second, with a well depth of -69.2 
cm~^ is located at ^ = 180° and R = 4.92 A. In Fig. 2b we show the minimum potential energy angular path, from 
which we see that a saddle point of -66.49 cm~^ is located at ^? = 53.5° and R = 4.29 A. 



B. Ground state calculations 



1. Variational Monte Carlo (VMC) method 

The ground state energies for the 0CS(p-H2) n clusters are first calculated using the VMC method. In this stochastic 
approach, one constructs a trial wavefunction 4'7'(R; {p}) which approximates the exact ground state wavefunction 
of the system, $o(R), for a given Hamiltonian H. Here, {p} denotes a set of adjustable parameters that control the 
shape of the trial wavefunction. The energy expectation value is then evaluated using a Monte Carlo integration. In 
most cases, \E't(R; {p}) ^o(R'), and, consequently, the VMC method provides an upper bound for the total energy, 

F < F({r.\) /^T(R;{p})g^T(R;{p})rfR 
The last equation can be rewritten in a convenient form for Monte Carlo integration as 

= /|*.(R;{p}) prfR 

1 " 

^-5^£;4R,;{p}) (4) 

k=l 

where n is the total number of configurations that are sampled from the probability density function | ^t(R; {p}) |^ 
using the Metropolis algorithm [16]. In eq. (4) i?j:,(Rfc; {p}) is the local energy 

ELiRk-, {p}) = *?'(R; {p})^*t(R; {p}). (5) 

In our implementation ^'t(R; {p}) is given by the usual generalized product form 

N N 

^T{-R;{p}) = l[^T{Rij;p')l[xT{Ri,^i;p) (6) 

where Ri, 'di and Rij are the same as in eq. (1). This type of trial function is symmetric with respect to p-&2 molecule 
exchanges and thus explicitly includes Bose statistics. 

For the (p-H2)-(p-H2) correlation, we used the two-parameter function 

^T(%)=exp(-|^-p^i?,,). (7) 

For the (p-H2)-0CS correlation, we started from an analytic fit to the ground state wavefunction of the A''=l (P-H2)- 
OCS complex calculated by use of the collocation method [17]. This function has the form 

XT{Ri,'&i) = exp{pirf= +p3[p9 +p4(cos7 -ps)'] Inr + [p5r'(cos7 - pg)' - Pio]e^^"^^''} (8) 

where r = \/y'^ + z^, y = Ri sini?j — pi2, z = Ri cos'&i — pu and C0S7 = z/r. 

Optimization of the parameters {p} for A'^ >1 is then performed by simultaneous minimization of the energy E{{p}) 
and its variance. This approach provides the best wavefunction \I/t(R; {p}) given a particular functional form [18]. 
The values of the optimized parameters are reported in Table I. 

In the present study, a VMC run typically consists of 5x10^ steps whose size is uniformly distributed and chosen 
such that the number of accepted moves is approximately mantained at half the number of the attempted ones. 
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2. Diffusion Monte Carlo (DMC) method 



Exact calculations of the ground state properties for a many-body quantum system can be made by use of the DMC 
method [19]. Wc employ here the rigid-body diffusion Monte Carlo (RB-DMC) scheme that is described in ref. [20]. 
In the following we provide only a summary of the main concepts and give details specific to the 0CS(p-H2) n system. 

Within the RB-DMC formulation we treat the OCS molecule as a rigid-rotor interacting with N spherical p-B.2 
molecules. The imaginary time (r = it/Ti) Schrodinger equation can then be written as: 

--[g--gre/]<&(R,r) (9) 

where Eref is the reference energy defining the zero of the absolute energy scale, and H is the Hamiltonian of the 

system 

In the r.h.s. of eq. (10) the first two terms are the translational and the rotational kinetic energy of the OCS molecule, 
respectively. Mq is the mass of the OCS and Bq its rotational constant, (jjx, 4'v describe the OCS rotations about the 
X- and y-axis in the molecule-fixed frame, where the ^;-axis is taken to lie along the OCS axis. The term — is 
the kinetic energy of the ith p-B.2 molecule and m is the p-B.2 mass. Finally, 1^(R) is the total interaction potential 
of eq. (1). 

It is well known that efficiency of the DMC method can be improved if one introduces a guiding function ^'t(R). For 

our 0CS(p-H2)jv system this leads to the importance-sampled RB-DMC equation for the product function /(R, r) = 
$(R, r)\E'T(R-), which differs from eq. (9) by the presence of additional drift terms: 

^2 N 2 N 

+^EV?/(R,-)-^Ev4/(R,r)F.(R)] 

i=l i=l 
-[EL{n) - Eref]f{Ii,T). (11) 



Here, £'i,(R) is the local energy of eq. (5), and 

f(*)(R) - Volnl vI/t(R) I' (12) 

Fi(R) = Vilnl *t(R) I' . (14) 

The last three equations define the quantum forces that direct the sampling of /(R, r) into regions where ^'t(R) is 
large. 

The steady state solution /(R, r oo) = $o(R)^t(R) of eq. (11) is obtained by use of a random walk technique. 

An ensemble of walkers {X^} with associated statistical weights {wk} is propagated in imaginary time from an initial 
arbitrary distribution using the short time approximation of the Green's function appropriate to eq. (11) [20]. Here, 
X is a vector in a (37V+5)-dimensional space. Two of these dimensions describe the OCS rotational motion and the 
remainder 3(A''+1) the translational motion of all the particles. The propagation from r to r + At is achieved by 
a combination of random gaussian displacements and systematic moves under the influence of the quantum forces. 
Detailed balance is ensured by use of a generalized Metropolis scheme at each time step [21]. 

The ground state energy Eq can then be calculated from averaging the local energy over the asymptotic distribution, 
according to 

^^'^ - //(R,r-oo)dR (^^^ 
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_ /$o(R)gvl/T(R)dR 

/$o(R)*T(R)dR ^ ^ 

= Eo, (17) 

where the last equaUty follows from acting with H to the left. 

Expectation values of coordinate operators A = ^(R) can also be obtained by averaging over /(R, r oo), but 
the integration leads here to a "mixed" expectation value 

,u I ^(R)/(R, T ^ oo)dR / <l>o(R)A(R)M/T(R)rfR 

'^^>'^''^- J f{K,T ^ oo)dR ~ /$o(R)*T(R)(iR ■ ^ ^ 

The "pure" expectation values can be retrieved after insertion of the factor [$o(R)/^t(R)] in eq. (18), i.e., 

,u _ /[^o(R)/^t(R)]A(R)/(R,t ^ oo)dR _ /<l>o(R)A(R)<l>o(R)rfR 

^ ^^""-^ /[$o(R)/*T(R)]/(R,T^oo)dR /$o(R)*o(R)ciR ^ ' 

In the present study the ratio [<I>o(R)/'I't(R)] is evaluated by use of the descendant weighting procedure as described 
in ref. [22]. All distributions shown here are computed with this method. We shall analyze specifically the two 
dimensional density distribution p{r, z) of the p-B.2 molecules around the OCS molecule. 
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^(^'^) = E(%r^^(^^-^))' (20) 

i 

where z = Rcosd is defined to lie along the OCS axis, and the polar radius r = Rsini} is perpendicular to this axis. 
We shall also show the pair angular distribution function: 

1 ^ 

Pi2{<t>) = ]v(ArrT) - '^))' (21) 

normalized such that 



/ 



Pi2(<^)# = l. (22) 



Here (jjij is the angle bc;twccn two P-H2 molecules and the OCS axis, in a plane perpendicular to the molecular axis. 
Therefore, P12 (^) is a three body distribution function integrated over all variables with the constraint of fixed (j). It 
describes the probability to find a IJ-H2 molecule at an angle cf) from another molecule that is arbitrarily located at 
(j)=0°. 

In all the DMC calculations presented here, the guiding function \E't(R) is given by eqs. (6-8). We use an ensemble of 
1000 walkers whose initial distribution is taken from a VMC run. A variable weight Wfe(r), with Wfc(O) = 1, is assigned 

to each walker. During the imaginary time propagation the weights can vary between WmAn and tUmax according to 
the mixed weights/branching scheme described in ref. [18], and the reference energy is updated continously according 
to the growth energy estimator 

p (^.j, lniy(r)-lniy(r-AT) 

Eref{T) = £^ref[T - At) — , [26) 

where W{t) = Ylk'^'ki'^)- Dependence on the size of the time step At has been carefully analyzed in order to 
eliminate any significant bias. A resulting optimal value of At=50 a.u. is used. After performing a thousand steps 
to equilibrate the initial ensemble, averages are calculated by performing block averaging on blocks consisting of 500 
steps, and subsequent averaging of these block averages. 



C. Excited state calculations 



1. Projection operator imaginary time spectral evolution (POITSE) method 

The POITSE method allows the calculation of excited state energies from a stochastic solution of the Schrodinger 
equation without imposing any nodal approximations. In this scheme, the excited state energies are extracted from 
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the two-sided inverse Laplace transform of an imaginary time correlation function k(t) that is obtained by a multi- 
dimensional Monte Carlo integration combined with zero temperature diffusion Monte Carlo sidewalks. The basic 
formulation of the method has been previously discussed [23] , and thus we present here only a brief summary of the 
main ideas, together with details specific to the present study. 
The analysis starts with the spectral function 

«('^) = I] I (V'o I i I V'n) I' S{Eo -En+u), (24) 

n 

where {| tpn)} and {En} are a complete set of eigenfunctions and eigenvalues for the Hamiltonian H, respectively. 
A is a local operator that projects, at least approximately, | ■^o) into the particular excited state of interest | Vn)- 
Performing a two-sided Laplace transform of eq. (24) results in the imaginary time correlation function 

p + co 

k{T) = / exp {—uJT)K{Lj)doj (25) 



= (V'o I iexp [-{H - Eo)t]A^ I Vo) (26) 
= 5]|(Vo|i|Vn)re-(^-^°)^ (27) 

n 

The matrix element of eq. (26) involves only the ground state and, therefore, it can be computed by using a multi- 
dimensional Monte Carlo (MC) integration over the distribution | V'o P and then performing a DMC sidewalk at each 
MC point. In most cases, however, the exact ground state | Vo) is not known. For practical purposes, one typically 
employs a trial function | '^t) and a reference energy Eref- These approximate | ipo) and Eq, respectively, as closely 
as possible. Use of a reference energy not equal to the exact ground state energy modifies the decay rate of all terms in 
eq. (27) by a constant factor Eref — Eq. It has been shown in ref. [23] that this dependence on the arbitrary reference 
energy can be eliminated by introducing the normalization factor 

(*T I exp [-{H - Eref)T] \ ^t)- (28) 

Thus, replacing | ^o), Eq in eq. (26) with | \1/t), Eref, respectively, and dividing by the additional normalization 
factor of eq. (28), leads to the desired time correlation function 

^ {^T \ Aexp[~{H- Eref)T]AU ^t) _ ^29) 

(*T I exp [~{H - Eref )T] I *t) 

From eq. (27) it is clear that the imaginary time decay of k(t) contains information about energy differences En — Eq. 
Therefore calculation of the numerical inverse Laplace transform of k(t) provides the corresponding spectral function 
k{ijj) whose peak positions correspond to the excitation energies En- 

The first step in the evaluation of ^(t) is a variational Monte Carlo walk in which an initial ensemble of 1000 
walkers distributed according to ^'j. is generated. The starting VMC ensemble is then propagated for 8000 steps of 
At=50 a.u. by a DMC sidewalk, during which eq. (29) is sampled. The trial function is given by eqs. (6-8), and the 
VMC and DMC implementations have been described in the previous sections. 

As well known, the numerical inversion of k(t) to obtain k{uj) is an ill-conditioned problem, especially when 
Monte Carlo noise is non- negligible and/or when the spectral function k{ui) contains multiple overlapping peaks of 
comparable intensity. Thus a judicious choice of the operator A is necessary to ensure that the time-dependence of 
K(r) is dominated by only one or few well-separated energy differences. 

In order to compute the rotational excitations of OCS inside clusters we employ here projectors A proportional 
to the molecular Wigner functions D-'^f.{a, (3,j). Here, a,/?, 7 are the three Euler angles that define the orientation 
of the molecular frame in the arbitrary space-fixed frame, and j, m, k are the OCS rotational quantmn numbers. We 
focus in particular on Dqq, with j ranging from 1 to 5. These projectors, except for a numerical constant, correspond 
to the rotational eigenfunctions of a linear molecule. They are also equivalent to the eigenfunctions of a symmetric 
top rotor with fc=0. They are functions only of the second Euler angle and are expressed in terms of the Legcndre 
polynomials, i.e., Dqq oc Pj(cos/3). Therefore, A accesses states in which the total angular momentum J is carried 
primarily by the OCS molecule, i.e. J ~ j. 

The inverse Laplace transform of k(t) is performed here using the implementation of the Maximum Entropy 
(MaxEnt) method as described in ref. [24]. Our use of this approach is identical to that employed in previous 
POITSE applications [23,25-27]. Since the maximum entropy analysis requires independent samples of k{t), the initial 



6 



configuration for each DMC sidewalk is taken from the VMC walk every 200 steps apart, to minimize correlations 
between successive sidewalks. In the present study, 8000-10000 independent decays arc required to produce a converged 
spectral function k{uj) for all the projectors defined before. We have also calculated the excitation energies by fitting 
the imaginary time correlation function with a sum of exponentials as was done in ref. [28], i.e., 

^(T)=^c„e-(^''-^°)^ (30) 

n 

Typically, three terms in this sum are sufficient to obtain convergent results. In eq. (30), the coefficient c„ is the 
spectral weight of the particular excitation En contributing to k(t). In order to discuss the convergence of our results 
we also analyse the quantity ej{T) defined as [28] 

e,(r) = -lln[K(^)(r)M")], (31) 

where c^^^ is the largest spectral weight obtained from the exponential fit of eq. (30). 

The analogies and the differences between MaxEnt and the exponential fit are discussed in more detail in section 
IIIB. 



2. Clamped coordinate quasiadiabatic diffusion Monte Carlo 

Excited state energies can also be calculated using the clamped coordinate quasiadiabatic diffusion Monte Carlo 
(ccQA-DMC) method proposed by Quack and Suhm [29]. We have previously applied this approach to analysis of 
small OCSC*Hc)Ar clusters [30] and the 0CS(p-H2) complex [31]. 

The key idea of this approximate approach is the assumption that the whole cluster follows the rotation of the 
dopant molecule. In the present implementation the instantaneous inertial tensor /(r) is computed at each step of 
the DMC propagation. Diagonalization of /(r) provides an effective centrifugal potential Vj_Ka.Kr(^)^ where J is 
the total angular momentum quantum number, and Ka, Kc are the pseudo-quantum numbers of an asymmetric rotor 
[32]. Adding Vj^Ka,Kc(^) to 1^(R) of eq. (1) leads to a modified Hamiltonian operator 

Hj^k^^kS^) = i?(R) + Vj.k^.kS^). (32) 

The RB-DMC method can then be used to solve the imaginary time Schrodinger equation for HjKa.K^i^), given 
a fixed set of J, Ka and Kc- Furthermore, because H(R) and Hj^Ka^K^T^) differ only by the centrifugal potential 
term, a correlated sampling scheme [33] can be used to calculate Eq together with all excitation energies Ej^k^,Kc in 
a single DMC run. This significantly reduces the statistical errors and also reduces the computational cost. 

Since the instantaneous inertial tensor is a coordinate operator, this is implicitly determined by the mixed distribu- 
tion "P-i '$0 (hcc cq.(18)). Consequently the ccQA-DMC excited states energies may be indirectly subject to some trial 
function bias. We discuss this and other possible sources of inaccuracy in the ccQA-DMC approach in Section IIIB 2. 

III. RESULTS 
A. Ground state 

Both the VMC and DMC ground state energies for the 0CS(p-Il2)Ar clusters are listed in Table II. We also report 
there the energy per J3-II2 molecule derived from the DMC results (column 4). The total energy decreases as a 
function of A'', indicating that the cluster is stabilized by adding P-H2 molecules. One can also see that up to A''=5 the 
difference between DMC and VMC energies is less than 5%. It subsequently increases for the largest sizes, illustrating 
the difficulty of designing accurate trial functions for clusters with > 6. The fact that A^=5 corresponds to some 
" magic" number is supported by inspection of the energy per single P-H2 molecule. This quantity first decreases with 
A^, reaching a minimum value at A'^=5 and then rapidly increases. A similar behavior was seen in the energetics for 
OCS(''He)Ar clusters [8]. 

The difference between clusters with N < 5 and those with A^ > 5 becomes more evident if one considers the 
para-hydrogen chemical potential /xjv = Eo{N) — Eo{N — 1). This is reported in Fig. 3. After decreasing up to 
^■=4, /ijv remains approximately constant at A^=4-5. It then markedly increases when A^=6. For N=7-8 /ijv is again 



7 



constant at a value slightly below the value found for A''=6. This indicates that the OCS(p-H2)6 cluster is relatively 

less stable than the others. 

The above energetic analysis suggests that some structural change might be occuring between A''=5 and N=6. This 
hypothesis is confirmed on inspection of Fig. 4 where we show 2-dimensional contour plots of the number density p{r, z) 
defined in eq. (20). Up to iV=5 the ground state structure of the clusters corresponds to a single ring of para- hydrogen 
molecules around the OCS axis. For these small sizes the P-H2 molecules are all located in the global minimum region 
of the 0CS-(p-H2) interaction (see Figs. 2a and 2b). From A'^=6 onwards, a second ring appears in the region close to 
the oxygen side. Note that this ring is not located in either of the two collinear local minima. Integration of p{r, z) 
for the A^=8 cluster shows that five P-H2 molecules arc still located in the first ring, while the other three are found 
in the second one. Note that while the first (lowest) local minimum of the 0CS-(p-H2) interaction potential is located 
at the sulfur end (i?=0°) and is about 20 cm~^ deeper than the second one at the oxygen end (see Fig. 2b), after the 
global minimum region is filled up with five P-H2 molecules, no P-H2 molecule is found at i? « 0°. This can be easily 
explained by considering the additional effect of the (p-H2)-(pH2) interaction. This provides an attractive potential 
that pulls the added P-H2 molecules closer to the first ring. Because of the saddle point at ■!?=53.5° the region close 
to the oxygen side {'&=0°) is preferred over the sulphur side (t?=180°). 

We also note that in the density plots, the peak corresponding to the first ring is separated from that of the second 
ring by a region where the density is essentially zero. This indicates that after the first ring is completed there is no 
delocalization of the other P-H2 molecules along the OCS axis. This is in contrast to what was previously found for 
^He in the analogous OCS('*He)]v clusters, where extensive helium delocalization is evident [8,30]. 

Although significant information about the cluster structures can be extracted from the number density p{r, z), this 
does not provide any information concerning the localization of the P-H2 molecules in planes perpendicular to the 
OCS axis. In the left panels of Fig. 5 we therefore show the pair angular distribution function defined in eq. (21), for 
clusters with N=2-6. In the right panels of the same figure we also report for comparison the corresponding results 
obtained for the same size OCS(^He)jv clusters. In this case we have used the He-OCS interaction potential of ref. 
[34] and the He-He potential proposed by Aziz and co-workers [35]. Pi2{4>) describes here the pair distribution of 
para-hydrogen (helium) in the ring located in the global minimum region of the 0CS-(p-H2) (OCS-He) interaction. 

At first sight it is evident that the OCS(p-H2)7v clusters show a stronger localization of the P-H2 molecules in this 
ring. It is also important to note that up to A'^=5 there is a zero probability to find two P-H2 molecules closer than 
their repulsive hard-core diameter (~ 3 A). This evidence together with the fact that this ring is finite and effectively 
one-dimensional implies that only cyclic permutations among P-H2 molecules are possible. For larger iV, instead, the 
probability to occupy sites outside the global potential well is finite (see Fig. 4), and this results in a finite value of 
Pi2{(t>) at 0=0°. 

The analysis of the pair angular distribution also provides some insight into the relative rigidity of the first ring 

as a function of the number of P-H2 molecules. For N=2, Pi2(0) shows two distinct peaks and thus the structure of 
the cluster can be reasonably described as a "rigid" P-H2 dimer that is located perpendicular to the OCS axis. When 
A''=3-4 the ring is instead "floppy", as evidenced by the broad central band in Pi2{4>)- In particular, we note that for 
N=4: there are not three distinct peaks corresponding to three p-B.2 molecules located at well-defined angle </> with 
the fourth at 0=0°. This means that the four P-H2 molecules are not arranged in a square configuration around the 
OCS, i.e., the ground state structure of OCS(p-H2)4 is not exactly symmetric. 

The most interesting aspect, however, is the structural transition that occurs when N=5. For this cluster size, in 
fact, Pi2(0) looks entirely different. This now shows four sharp peaks, corresponding to four well-localized P-H2 with 
the reference P-H2 at 0=0°. The OCS(p-H2)5 structure can thus be rationalized as a "rigid" ring of para-hydrogen 
molecules around the OCS axis. In contrast, the OCS(^He)jv clusters appear to be delocalized in the degree of 
freedom and, consequently, very floppy for all the sizes shown in Fig. 5. 

B. Excited states 

1. POITSE results 

In Section III. A we have seen that the first few para-hydrogen molecules form a ring around the OCS axis. However, 
as discussed there, for < 5 the ring is not complete and the P-H2 are not symmetrically arranged (see Fig. 5). This 
implies that the molecular axis is not the axis of symmetry of the whole cluster which, as a result, has an asymmetric 
structure. 

As well known, the rotational eigenfunctions for an asymmetric rotor are given by a linear combination of the 
Wigner functions [36]. In POITSE calculations, however, we use projection operators that are proportional only to 
one of these functions, namely Z)QQ(a, /J, 7). Furthermore, the three Euler angles are defined with respect to a frame 
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fixed on the OCS molecule that, as noted above, does not coincide with the principal axis frame of the whole cluster. 
As a consequence of this, the excited state wavefunction obtained after the action of Z?QQ(a, /3, 7) on the ground state 
wavefunction, is not an eigenstate of the Hamiltonian. Therefore, it can have a non zero overlap with excited states 
different from those we are considering here. In a symmetric top limit, these states correspond to states with J=j 
and K=k=0, where J is the total angular momentum quantum number, and K, k are the quantum numbers for the 
total angular momentum projection along the space-fixed and the body-fixed z axis, respectively. 

That this is indeed the case we are facing here becomes clear if one looks carefully at Fig. 6. In the right panels, 
we show the spectral function obtained from the Maximun Entropy algorithm for N=l-A and J=l-3. In the left 
panels we show the corresponding quantity ej(T) (solid lines) of eq. (31), and the excited state energies (dotted lines) 
obtained from the exponential fits, eq. (30). 

It is evident that for J=l the quantity ej(r) is constant in time for all the four clusters, indicating that the projection 
operator Dqq accesses primarily one eigenstate of the Hamiltonian. However, for J=2 and J=3, (-,j{t) first increases, 
reaches a maximum value corresponding to the excitation with the largest spectral weight, and then decreases, as a 
function of r. Furthermore, the larger is J, the more pronounced is this last decrease. This is the consequence of 
the fact that, because of the asymmetric structure of these smallest clusters, the projection operators Dqq can access 
different rotational excited states. Thus the corresponding imaginary time correlation function contains more than 
one eigenergy and €{t) varies as these eigenenergies are accessed. It is also important to note that when A^=4, the 
OCS axis more closely approximates the axis of symmetry for the whole cluster, and hence e(T) becomes constant in 
time also for the higher values of J. 

From these considerations, it follows that the inversion of the correlation function k(t) is technically more difficult 
for the highest values of J and for N < 3. In fact, while for J=l the result of inversion is independent of the time r, 
for J=2 and J=3 MaxEnt provides a stable spectral function with a narrow peak only when the inversion range is 
restricted to the time interval before e(T) decreases. The same sensitivity to time interval is found when the inversion 
is perfomcd using the exponential fit. Identical results are obtained from the two inversion procedures. However, in 
the case of the exponential fit, the optimal fit results in a very large value (x^ Ri 5 x 10^), indicating that caution 
must be used when dealing with this procedure. 

From N=5 onwards, we have seen that a rigid ring of five P-H2 molecules is formed around the OCS axis and thus 
the molecular axis can be reasonably identified with the axis of symmetry of the whole cluster. As a consequence of 
this, the I?oo's become "good" projection operators each of which accesses only one eigenstate of the Hamiltonian. 
The inversion of k,{t) is thus simpler for these sizes and we are able to extract the excited state energies for J ranging 
from 1 to 5. In Fig. 7 we report, as an example, e,7(r) and k{ll!) for A^=5 and N=8. From the use of the Maximum 
Entropy method we have found that, while the peak positions are independent of the length of the time interval, 
narrower peaks can be obtained if only the central part (from r=1.0xlO^ a.u. to r=3.0xl0^ a.u.) of the time decay 
is considered. For these cluster sizes the exponential fit also provides identical results as MaxEnt, but again the 
associated x^ is very large (see above). 

This detailed analysis of the imaginary time dependence of €{t) shows that the complexes become significantly 
symmetric for N > 5. 

2. ccQA-DMC results 

In Fig. 8 we show the rotational excitation energies Ej^Ka.,Ka with J=l-3 calculated by use; of the c;cQA-DMC 
approach for A^==l-8. Each Ej^Ka,Kc exhibits a monotonic decrease as N increases. The patterns resemble those of 
slightly asymmetric prolate rotors. Prom A''=4 onwards, the energy differences between those excited states that are 
degenerate in the limit of a perfect symmetric prolate rotor become smaller and smaller. This indicates that, as the 
POITSE results have already shown (Figs. 6 and 7), the clusters become more symmetric. Note that A^=5 is the 
least prolate structure and most symmetric structure in this size range. This can be understood as a consequence of 
the fact that this cluster possesses a single complete ring. It is also possible to see from Fig. 8 that the structural 
transition between A^=5 and N=6 evident in Fig. 4 manifests itself in a change of the slope in Ej^Ka,,K^- This is 
particularly evident for the lowest excited states (bottom panel). 

In Fig. 9 we compare the exact POITSE with the approximate ccQA-DMC energy levels for J=l-3. Because 
the cluster structures arc those of slightly asymmetric prolate rotors, the excitations extracted from the molecular 
projection operators in POITSE correspond to those of states with J=j, i^a=0 and Kc=l- The two sets of results 
show generally similar behavior but there are nevertheless several significant differences. For ^"=1 and N=2 the 
POITSE values are systematically higher than the corresponding ccQA-DMC values, indicating that for these sizes 
the P-H2 molecules are not as rigidly bound as assumed by ccQA-DMC. For A^=3-4, in contrast, the two sets of 
results are identical for J=l but differ increasingly as J increases, with the greatest discrepancies being seen at J=3. 
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Such a pattern of differences is identical to what would be expected for excitations lowered by centrifugal distortion 
below the corresponding rigid molecule values. This can be rationalized by considering that for these two sizes the 
pair angular distribution function has shown that the structures are quite floppy (Fig. 5). Consequently, for higher J 
values the p-B.2 density can be increasingly easily distorted. The rotational energies indicate that iV=3 is therefore the 
floppiest complex. This effect consistently disappears for A''=5 and A''=6 when a complete ring of five molecules 
is formed around the OCS axis. At these sizes the POITSE and ccQA-DMC energy levels are essentially coincident, 
implying that the true energy levels (POITSE) correspond now to a rigid complex in which all P-H2 molecules are 
rigidly coupled to the OCS rotation (ccQA-DMC). We can rationalize these different behaviors for the pairs A''=3-4 
and A^=5-6 by a coupling of the ring "breathing" mode to the overall rotation in clusters with an incomplete ]3-H2 ring 
at the global potential minimum (see Fig. 4). Such coupling results in a centrifugal distortion of the cluster, due to 
expansion of this floppy para- hydrogen ring. The vibrational "breathing" mode is relatively facile for an incomplete 
ring. When this is filled with 5 pH2, the (p-H2)-(p-H2) interactions effectively lock the para-hydrogen molecules into 
a more rigid structure (see Fig. 5 and effectively eliminate centrifugal distortion due to ring "breathing". 

For N—7-8 quite different behavior is seen. Here the POITSE values lie slightly below the ccQA-DMC energies for 
all J values, including J=l. This does not correspond to the pattern of a centrifugal distortion away from a rigid 
complex. The systematically lower POITSE energies for these cluster sizes are somewhat surprising. In fact, because 
of the rigid coupling assumptions implicit in the approach, the ccQA-DMC energies are expected to provide a lower 
limit, at least for the lowest value of J where non-rigidity effects are less important. Resolution of these differences 
lies in the implicit dependence of ccQA-DMC on a "mixed" ground state density. Fig. 10 compares the dependence 
of the "mixed", eq. (18), and "pure", eq. (19), ground state densities for 7V=l-8. We sec that while for A^=l-6 these 
densities are very similar, for N=7-8 the "mixed" state density overestimates the true density in the region between 
the two rings, at the expense of the second ring density. This results in an underestimate of the moments of inertia 
Ib and Ic, and hence in an overestimate of the rotational energies obtained from ccQA-DMC. 

3. Rotational constants 

The rotational excited state energies obtained from the POITSE method are fitted using a symmetric top Hamil- 
tonian. For states with K=0, this becomes 

E{J) = BavgJiJ+l)-Dj\j + lf. (33) 

Because the ground state cluster structures correspond to those of slightly asymmetric prolate top rotors, at least 
for the smallest sizes (Figs. 4 and 5), the averaged rotational constant Bavg refers to {B + C)/2. Both Bavg and the 
distortion constant D are shown as a function of N in Fig. 11. 

The averaged rotational constant Bavg decreases as a function of N, showing a marked change in its slope between 
A''=5 and A'"=6, exactly when the P-H2 molecules begin to solvate the oxygen side of the OCS. As we have discussed 
in the previous section, at A''=6, one P-H2 molecule occupies a site closer to the OCS axis and more distant from the 
perpendicular rotation axes. Consequently the moment of inertia increases providing a larger reduction of Bavg- It is 
also important to note that a similar magnitude decrease of Bavg from A^=l to N=3 and from A'^=6 to A'^=8 is found. 
This suggests that the second ring around the OCS axis has similar azimuthal distribution features as the first one. 

The distortion constant _D is ^ 4 x 10~^ cm^^ for A^=l, appearing somewhat larger than that for the OCS-'*He 
complex [5], although the large error bars make a meaningful comparison impossible. D is smaller for N=2, nearly 
zero, providing more evidence that the first two p-B.2 molecules constitute a near "rigid" dimer interacting with the 
OCS. Beyond N=2, the distortion constant shows a maximum value at A''=3 that, as we have already discussed above 
(Section III.B.2), corresponds to the floppiest structure. This is followed by a slightly lower, but still significant, value 
at A^=4 and a subsequent decrease with N, because of the increasingly rigid pac;king of P-H2 molecules around the 
OCS. The relatively substantial values of D for A''=3-4 are due to the distortion of the ring deriving from coupling 
to the ring "breathing" mode, described in Section III.B.2 above. As noted there, for A''=5-6 five hydrogen molecules 
are now locked in the ring, which will push this breathing mode higher in energy and reduce its amplitude, thereby 
reducing the effects of its coupling to the cluster rotation. The continued decrease of D for larger N {N=7-8), is very 
interesting in this context. For N=7-8 the second shell contains 2 and 3 para-hydrogen molecules respectively. One 
might therefore expect that for the A^=8 chister, the "breathing" mode of the second 3-moleculc ring would similarly 
couple to the rotation and cause some centrifugal distortion. However, interaction of the 3 P-H2 molecules in this 
secondary ring with the 5 molecules in the primary ring stabilizes and localizes the secondary ring, reducing the effects 
of any rotational coupling to its breathing mode. Consequently, there are no appreciable centrifugal distortion effects 
due to growth of the second ring. 
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As a result of the approximations made in the method, the ccQA-DMC calculations can provide estimates for all 
three rotational constants {A, B, C) within its rigid coupling assumption. These are reported in Table III together 
with a comparison of the resulting values of ^avg with those obtained from POITSE. The ccQA-DMC values confirm 
the increasingly symmetric prolate structures as N increases. For A''=5 the rotational constants B and C differ by 
only 0.013 cm~^, while for A'' > 6 this difference is less than 0.01 cm~^. We interpret these as essentially symmetric 
top structures. The differences between the two sets of results for ^avg arc very small for all cluster sizes studied here, 
indicating that the rigid coupling approximation made in the ccQA-DMC method is generally a good approximation. 
Moreover, the differences seen here are smaller than those found in the analogous OCS(^He) m clusters [10,28] , implying 
that the complexes with para-hydrogen are generally more rigid than those with helium. However, the ccQA-DMC 
approach does rely on several approximations that can cause small inaccuracies. We have already noted one problem, 
namely the implicit dependence on mixed densities that affected the ccQA-DMC energies for the N=7-8 clusters. 
Another general problem is that the ccQA-DMC method overestimates the degree of asymmetry of the clusters, since 
the averaging of fluctuations can dominate the results even for highly symmetric structures. Thus for the symmetric 
OCS(p-H2)5 cluster, for which the ground state structure is clearly symmetric in the azimuthal degree of freedom 
(Fig. 5), ccQA-DMC nevertheless yields B ^ C. This suggests that some small quantitative inaccuracies might be 
present for A/" > 5. 

IV. RELATION TO EXPERIMENTAL MEASUREMENTS 

In the following, the present results will be discussed in light of the existing experimental studies of 0CS(p-H2)jv 
complexes. 

For N—1 a direct comparison with the data reported in rcf. [5] is possible. The high resolution spectrum for the 
0CS(p-H2) complex was observed and analysed using the conventional asymmetric rotor Hamiltonian in the a-reduced 
form of Watson [37]. From that analysis all the rotational parameters were obtained. Structural properties as the 
intermolecular distance R and the Jacobi angle -d were also derived. As discussed before, the POITSE calculations 
provide rotational excited energies corresponding only to states with J=j, AT^^O and Kc=J ■ Therefore, from this set 
of calculations it is not possible to extract all the rotational constants. In Table IV we report the energies calculated 
from POITSE and make a comparison with the corresponding experimental values. Excellent agreement is found for 
all the three excited states considered here. The values for {R) and {§) calculated by DMC and reported in Table V 
also show a very good agreement with the corresponding experimentally determined quantities. A complete analysis 
of the ccQA-DMC results for A''=l was previously reported [31]. 

For A' > 1 no experimental results are available for the pure 0CS(p-H2)Ar clusters. Infrared spectra for these 
clusters were however observed in pure "^Hc and in mixed ^Hc/'^He droplets [1-3]. Although a direct comparison with 
those results is not possible because it is difficult to estimate the effect of the surrounding helium droplet, we may 
still make some qualitative comparisons. The DMC calculations (ground state and ccQA-DMC) provide asymmetric 
structures for clusters with N < 5 (Figs. 5-8, and Table III). This is confirmed by analysis of POITSE results for 
the rotational excited state energies (Figs. 6-7 and discussion in Section III.B.l). The present results are thus in 
good agreement with the experimental observations that allow a fit of the rotational spectrum using a symmetric top 
Hamiltonian only for A^ > 5 [2]. For A^ = 5, the rigidity of the single ring structure seen in the DMC structures and 
in the POITSE energies provides justification for the permutation symmetry arguments made to explain the lack of 
a Q-branch in the spectrum at this size [2]. 

It is important to point out that the experimental data in helium droplets suggest that the principal ring located 
in the minimum region of the 0CS-(p-H2) interaction should contain six P-H2 molecules [3]. This appears to be 
in contradiction with the present results which show the first ring containing only 5 para-hydrogen molecules. The 
difference can be easily explained when the effect of the helium density in the droplet is taken into account. It has 
recently been demonstrated with path integral Monte Carlo calculations that in the mixed cluster 0CS(p-H2)M(^He)jv, 
six j)-H2 molecules are effectively arranged in the primary ring [38] . A ring of six P-H2 molecules is also foimd in the 
pure 0CS{p-il2) N clusters with A^ >8 [10]. The number of para-hydrogen molecules in the primary ring depends on 
the balance of (p-H2)-(p-H2), (p-H2)-0CS, and (p-H2)-^He interactions. 

V. SUMMARY AND CONCLUSIONS 

The work reported here presents a detailed study of the energetics and structures of the OCS(p-H2)7v clusters with 
A^=l-8. The rigid-body diffusion Monte Carlo method has been employed to calculate the ground state properties. 
Analysis of the total energy as a function of the cluster size N shows that the clusters are stabilized by adding P-H2 
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molecules. Analysis of the chemical potential indicates that clusters with A^=4-5 arc relatively the most stable. 
For these two cluster sizes the calculation of the P-H2 density distribution shows that a single ring of P-H2 is formed 
around the OCS axis and is located in the global minimum region of the 0CS-(pH2) interaction. 

Inspection of the p-B.2 angular distribution inside this ring provides important informations about its rigidity. For 
N=2 the cluster is effectively constituted by a near "rigid" P-H2 dimer bound to the OCS molecule. The clusters with 
N=3-4: appear to be floppy. When A^=5 the ring is complete and the P-H2 molecules arc again "rigidly" arranged 
around the molecular axis. Consequently, the OCS(]>H2)5 cluster corresponds to a symmetric structure. 

The rotational energies for excited states with J w j have been also computed as a function of the cluster size N by 
using the POITSE method. Comparison with results obtained from a rigid coupling approximation, the ccQA-DMC 
approach, shows that over this size range the clusters of OCS with J3-H2 are to a first approximation fairly "rigid", 
showing less deviations from near "rigid" behavior for N=6-8 than do complexes of OCS with helium [10,28]. It has 
been shown that the detailed variations in rotational spectra of the clusters with A''=l-8 can be readily explained in 
terms of their structural features. Thus, for < 4 the rotational spectra correspond to those of slightly asymmetric 
prolate rotors with significant centrifugal distortion for A^==3-4 that is assigned to coupling of the rotational motion 
to a "breathing" mode vibration of the partially filled para-hydrogen ring. This coupling is suppressed at A^=5 when 
the ring is complete and the (j>-H2)-(p-H2) interactions effectively lock the para- hydrogen molecules into a more rigid 
structure. For N > 5 prolate symmetric top spectra are therefore obtained with negligible distortion constants. Some 
quantitative inaccuracies in the ccQA-DMC method were seen to arise when the second ring is growing (7V=7-8) and 
when the cluster is highly symmetric (A^=5). 

Comparison with the experimental data for the A=l 0CS(p-H2) complex shows excellent quantitative agreement 
for both the energy levels and the structural parameters. For A^ > 2 no experimental data are available for the pure 
0CS(p-H2)Ar clusters. A qualitative comparison of the main features of the cluster structures with conclusions drawn 
from the experimental results for the 0CS(p-H2) n clusters in the mixed '*Hc/"^Hc droplets shows good agreement, with 
the main difference being the helium-induced modification of the A'^=6 structure from a 5-membered to 6-membered 
ring. The rigidity of the highly symmetric A''=5 structure provides strong support for the arguments that excited 
states of this are restricted by nuclear permutation symmetry, resulting in a lack of Q-branch in the OCS rotational 
spectrum at low temperatures [3] . 
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TABLE I. Variationally optimized parameters for the trial functions and xt of eqs. (6-8). All values in a.u. 



N 


Pi 


P2 


P3 


Pi 


P^> 




Pr 


Ps 


Pfl 


Pio 


Pii 


Pl2 


P'l 


p'2 


1 


-4.35 


0.793 


6.12 


-0.22 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.625 


0.13 






2 


-4.33 


0.793 


6.12 


-0.22 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.625 


0.13 


6550 


0.05 


3 


-4.33 


0.793 


6.12 


-0.22 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.625 


0.13 


6550 


0.05 


4 


-4.33 


0.793 


6.12 


-0.22 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.625 


0.13 


6550 


0.05 


5 


-4.33 


0.793 


6.12 


-0.22 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.625 


0.13 


6550 


0.05 


6 


-4.33 


0.793 


6.12 


-0.16 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.725 


0.13 


6550 


0.05 


7 


-4.33 


0.783 


6.12 


-0.05 


-0.16 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.800 


0.13 


6550 


0.05 


8 


-4.33 


0.773 


6.12 


0.10 


-0.20 


6.38 


0.958 


-0.188 


0.94 


1.095 


-0.800 


0.13 


6550 


0.05 



TABLE IL VMC and DMC energies for OCS(j9-H2)jv clusters with N=l-8. All values in cm"^ 



?,7 ,.^'^J(■ i^DMc ,•D^J('l^r 
f2 £1" £io ^0 

1 -74.765±0.003 -74.927±0.002 -74.927±0.002 

2 -152.088±0.007 -154.206±0.007 -77.103±0.003 

3 -231.89±0.01 -235.09±0.01 -78.363±0.003 

4 -312.98±0.02 -319.50±0.02 -79.875±0.004 

5 -385.1±0.03 -403.57±0.03 -80.715±0.005 

6 -427.1±0.04 -460.4±0.05 -76.740±0.008 

7 -473.1±0.05 -521.7±0.06 -74.527±0.009 

8 -5;!i.7±l).()5 -r)8:!.5=().()8 -72.91=0.01 



TABLE IIL Rotational constants for OCS(p-H2)jv clusters with iV=l-8. Bavg is defined as {B + C)/2. All values in cm 
The numbers in parentheses indicate the statistical error in units of the last digit. 





ccQA-DMC 


POITSE 


N 


A 


B 


C 


^avg 


^avg 


1 


0.7582(1) 


0.1902(2) 


0.1508(2) 


0.1705(2) 


0.179(4) 


2 


0.36756(5) 


0.1601(1) 


0.1349(1) 


0.1475(1) 


0.152(5) 


3 


0.23735(3) 


0.14445(7) 


0.11680(6) 


0.13062(6) 


0.136(2) 


4 


0.17477(3) 


0.12556(6) 


0.10699(5) 


0.11627(6) 


0.120(3) 


5 


0.13761(2) 


0.11029(5) 


0.09797(4) 


0.10413(5) 


0.108(1) 


6 


0.12533(3) 


0.08984(7) 


0.08052(5) 


0.08518(6) 


0.088(1) 


7 


0.11050(3) 


0.07631(7) 


0.06806(5) 


0.07218(6) 


0.068(1) 


8 


0.09652(3) 


0.06556(8) 


0.05923(5) 


0.06240(7) 


0.0566(9) 



TABLE IV. Comparison between POITSE and experimental results [5] for the 0CS(p-Il2) rotational excited state energies 
Ej^Ka,Ko with J=j, Ka=0 and Kc=J- All values in cm"'^. The numbers in parentheses indicate the statistical error in units 
of the last digit. 



J 




Ko 


POITSE 


Exp. 


1 





1 


0.355(6) 


0.3534 


2 





2 


1.06(2) 


1.0574 


3 





3 


2.08(8) 


2.1068 
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TABLE V. Comparison between DMC and experimental results [5] for the intermolecular distance (R) and the Jacobi angle 
{'&) (see Fig. 1). {R) is in A and {•&) in degrees. The numbers in parentheses indicate root mean square deviations. 





DMC 


Exp. 


(R) 


3.794 (0.365) 


3.719 




105.7 (9.2) 


110.8 
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FIG. 2. (a) Contour plot of the 2-dimcrisional OCS-(p-H2) potential energy surface obtained by averaging the 4-dimensioual 
potential of Higgins et al. [15] over the p-H2 angular variables {§' R and are the Jacobi coordinates defined in Fig. 1. 
The inner isoline around the global minimum at Ji=3.35 Aand i?=105° correponds to -140 cm~^. The spacing between two 
consequent isolines is 20 cm~^. (b) Minimum potential energy path in the angular coordinate i?. 
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FIG. 5. Pair angular distribution function Pi2(0), eq. (21), for N=2-6. Left panels: 0CS(pH2)Ar clusters, right panels: 
OCSC'He)jv clusters. 
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FIG. 6. e(r) calculated by means of exponential fit eqs. (30,31) and spectral function k{uj) eq. (24) calculated using Maximum 
Entropy method for N=l-4. 
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FIG. 7. e(r) calculated by means of exponential fit eqs. (30,31) and spectral function k{uj) eq. (24) calculated using Maximum 
Entropy method for A'^=5 and A'^=8. 
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FIG. 10. Angular profiles of the P-H2 densities for the OCS(p-H2)jv clusters with Af=l-8. Dotted lines: "mixed" densities 
from eq. (18), solid lines: "pure" densities from eq. (19). For A''=l-5 these densities are essentially coincident. 
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FIG. 11. Rotational constant Bavg (bottom panel) and distortion constant D (top panel) obtained from POITSE calculations 
for OCS(p-H2)jv clusters with N=l-8. 
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